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Light quark masses from unquenched lattice QCD 
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We calculate the light meson spectrum and the light quark masses by lattice QCD simulation, 
treating all light quarks dynamically and employing the Iwasaki gluon action and the nonperturba- 
tively 0(a)-improved Wilson quark action. The calculations are made at the squared lattice spacings 
at an equal distance a 2 ~ 0.005, 0.01 and 0.015 fm 2 , and the continuum limit is taken assuming an 
0(a 2 ) discretization error. The light meson spectrum is consistent with experiment. The up, down 
and strange quark masses in the MS scheme at 2 GeV are m = (m u + md)/2 = 3.55lg'2g MeV and 
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_ 6 { MeV where the error includes statistical and all systematic errors added in quadra- 



m s = 90. i: 

ture. These values contain the previous estimates obtained with the dynamical u and d quarks 
within the error. 
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The masses of light quarks are fundamental parame- 
ters of QCD. They cannot be measured experimentally 
since quarks are confined in hadrons. Lattice QCD en- 
ables calculations of hadron masses as functions of quark 
masses, and hence allows a determination of the quark 
masses from the experimental hadron masses. This ap- 
proach has been successfully applied, first in quenched 
QCD [1] and then in Nf = 2 QCD where degenerate 
up (u) and down (d) quarks are treated dynamically 
These studies have revealed that the light quark mass val- 
ues are significantly reduced by dynamical u and d quark 
effects. In this article, we present our attempt to deter- 
mine the quark masses in Nf — 2 + 1 QCD where the 
heavier strange (s) quark is also treated dynamically. We 
wish to examine to what extent the dynamical s quark 
affects the light quark masses. We determine the quark 
masses in the continuum limit and estimate all possible 
systematic errors. We also calculate the prerequisite light 
meson spectrum. A similar attempt has been made by 
the MILC Collaboration [|. 

We adopt the Iwasaki RG gauge action 0] and the 
clover quark action with the improvement coefficient csw 
determined nonperturbatively for the RG action [Bj] . The 
choice of the gauge action is made to avoid a first-order 
phase transition (lattice artifact) observed for the plaque- 
tte gauge action [6( . We employed the Wilson quark for- 
malism because we prefer an unambiguous quark-flavor 
interpretation over the computational ease of the stag- 
gered formalism adopted by the MILC collaboration [7| ■ 

Configurations are generated at three values of the cou- 
pling (3 = 6/g 2 — 2.05, 1.90 and 1.83 corresponding to the 



squared lattice spacing a 2 ~ 0.005,0.01 and 0.015 fm 2 , 
with the physical volume fixed to about (2.0fm) 3 . At 
each (3, we perform simulations for 10 quark mass com- 
binations using a combined algorithm [8] of the Hybrid 
Monte Carlo (HMC) for the degenerate u and d quarks 
and the polynomial Hybrid Monte Cairo (PHMC) for the 
s quark. Table U summarizes the simulation parameters. 

The meson and quark masses at the simulation 
points are determined from single exponential correlated 
X 2 fits to the correlators (P(i)P(O)), (V(t)V(O)) and 
(Ai(£)P(0)}, where P, V and denote pseudoscalar, 
vector and nonperturbatively (3(a)-improved axial- 
vector current operators, respectively. We use an ex- 
ponentially smeared source and a point sink, and mea- 
surements are made at every 10 HMC trajectories in the 
Coulomb gauge. For the pseudoscalar sector, (P(i)P(O)) 
and (Ai(i)P(O)) are fitted simultaneously ignoring cor- 
relations among them. Errors are estimated by the jack- 
knife method with a bin size of 100 HMC trajectories; 
errors do not increase for larger bin sizes. 

Chiral fits are made to the light-light (LL), light- 
strange (LS) and strange-strange (SS) meson masses 
simultaneously ignoring their correlations, using a 
quadratic polynomial function of the sea quark masses 
{m u ,md,m s ) and valence quark masses (m V aii, m va i2) in 
mesons; 

/(M S ,M V ) (1) 
= A + P s trM s + BytrM v + D sl /trM s trM v 

+C 5 itrM s 2 + C S2 (trM s ) 2 + CyitrM 2 + Cy 2 (trM v ) 2 , 

where / = m 2 PS or vector meson mass my, Mg = 
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TABLE I: Simulation parameters; L 3 x T is the lattice size, {n u d, k s ) is the hopping parameter combination, 1/St is the number 
of molecular dynamics steps in one trajectory, jV po iy is the PHMC polynomial order, and traj. is analyzed trajectory length. 
Pseudoscalar vector mass ratios ^ES- are also listed for light-light (LL) and strange-strange (SS) mass combinations. 



P = 1.83, L [i x T = 16 3 x 32, c sw = 1-761 



St N poly traj. ^(LL) ^(SS) 


St iV P oi y traj. ^f(LL) ^f(SS) 


0.13655 0.13710 1/80 80 7000 0.7772(13) 0.7522(15) 
0.13710 1/85 80 7000 0.7524(21) 0.7524(21) 
0.13760 1/100 100 7000 0.7076(18) 0.7414(17) 
0.13800 1/120 110 8000 0.6629(22) 0.7365(16) 
0.13825 1/140 120 8000 0.6213(24) 0.7343(15) 


0.13655 0.13760 1/90 110 7000 0.7769(14) 0.7235(19) 
0.13710 1/100 110 8600 0.7448(14) 0.7128(16) 
0.13760 1/110 120 8000 0.7033(18) 0.7033(18) 
0.13800 1/120 130 8100 0.6525(23) 0.6941(20) 
0.13825 1/150 150 8100 0.6083(32) 0.6884(21) 


P = 1.90, L A x T = 20 J x 40, csw = 1-715 


Sr N poly traj. ^f(LL) ^(SS) 


Sr N poly traj. ^s(LL) ^s(SS) 


0.13580 0.13580 1/125 110 5000 0.7673(15) 0.7673(15) 
0.13610 1/125 110 6000 0.7435(18) 0.7647(17) 
0.13640 1/140 110 7600 0.7204(19) 0.7687(15) 
0.13680 1/160 110 8000 0.6701(27) 0.7673(17) 
0.13700 1/180 110 7900 0.6390(22) 0.7691(15) 


0.13580 0.13640 1/125 140 5200 0.7667(16) 0.7211(21) 
0.13610 1/125 140 8000 0.7444(15) 0.7182(17) 
0.13640 1/140 140 9000 0.7145(16) 0.7145(16) 
0.13680 1/160 140 9200 0.6630(21) 0.7127(17) 
0.13700 1/180 140 7900 0.6243(28) 0.7102(20) 


P = 2.05, L A x T = 28 a x 56, c sw = 1.628 


Sr N poly traj. ^f(LL) ^f(SS) 


«. St iV P oi y traj. ^f(LL) =£s(SS) 


0.13470 0.13510 1/175 200 6000 0.7757(26) 0.7273(29) 
0.13510 1/195 200 6000 0.7316(24) 0.7316(24) 
0.13540 1/225 200 6000 0.6874(30) 0.7395(23) 
0.13550 1/235 200 6500 0.6611(34) 0.7361(25) 
0.13560 1/250 200 6500 0.6337(38) 0.7377(28) 


0.13470 0.13540 1/175 250 6000 0.7790(23) 0.6821(32) 
0.13510 1/195 250 6000 0.7341(29) 0.6820(39) 
0.13540 1/225 250 6000 0.6899(34) 0.6899(34) 
0.13550 1/235 250 6500 0.6679(45) 0.6899(43) 
0.13560 1/250 250 6500 0.6361(47) 0.6852(46) 



p=1 .90, L xT=20 x40, polynomial 



|3=1 .90, L xT=20 x40, polynomial 
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TABLE II: Lattice spacings in fm units. 
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FIG. 1: Chiral fits of meson masses with m^ WI at P — 1.90. 

diag(m u , m d , m s ), M v = diag(rn va ii, rn va .\ 2 ), and "tr" 
means the trace of matrices. In the fits, we use 
the axial-vector Ward identity quark mass m q = 
lim t ^oo(04^4(«)P(O))/(2<P(i)P(O)>) and set A = B s = 
Csx — Cs2 = for rrip S . These fits reproduce mea- 
sured data well, as illustrated in FigfTJ with reasonable 
X 2 /d.o.f. of at most 1.36. 

The physical quark mass point and the lattice spac- 
ing are determined from the experimental values of ir°, 
p° and K (if-input) or 7r°, p° and tf> (0-input) meson 
masses. Taking the p° mass as input may cause a large 
systematic error, because the p — > 7T7r decay mode is 
not open for our mass range (the lightest pion mass in 
this simulation ~ 620 MeV) and hence chiral extrapo- 
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lation of my for lighter quarks may be quite different 
from our fits. In order to estimate this uncertainty, we 
also check another combination [ir°,K, <£>]. We assume 
the ideal mixing for the vector isosinglets. Since our sim- 
ulation is made with degenerate u and d quarks, we con- 
sider the isospin averages m,j> — {{m 2 K± + m 2 K0 ) /2} 1 / 2 



and m^, — [m K *± + m K *o)/2 and predict the aver- 
age light quark mass m = (m u + m^)j1. The electro- 
magnetic (EM) effects, not included in our simulations, 
are removed from the m K ± above using Dashen's the- 
orem (m 2 K± - m 2 K0 ) E M = {ml ± - m 2 ) E xp- The 
isospin breaking effects and the EM effects for other 
mesons we consider are expected to be small and thus are 
not considered. The experimental values we use are taken 
from the PDG booklet [U; m,o = 0.1350GeV, m n ± = 
0.1396GeV, m K o = 0.4976GeV, m K ± = 0.4937GeV, 
m p o = 0.7755GeV, m K ,o = 0.8960GeV, m K ,± = 
0.8917GeV and m = 1.0195GeV. Lattice spacings (Ta- 
ble[Tl| for the K- and <jy inputs are consistent, while those 
for the [it, K, 0]-input are slightly smaller by at most 7%. 

An agreement of the meson spectrum with experiment 
is a necessary condition for a reliable estimate of the 
quark masses. To confirm this, we extrapolate the me- 



3 



TABLE III: Meson masses in the continuum limit (in MeV 
units), compared to experiment. The EM effect is subtracted 
using Dashen's theorem. 
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FIG. 2: Continuum extrapolation of meson masses for Nf = 
2 + 1 QCD (circles) , compared to experiment (stars) and re- 
sults in Nf = 2 (squares) and Nf = (triangles) QCD [2]. 



son masses linearly in a 2 , because our action is 0(a) 
improved and data are well fitted, as shown in Figf^l 
with small x 2 /d.o.f < 1.4. The masses in the continuum 
limit, summarized in Table HTT1 are consistent with exper- 
iment with at most 2.9a deviation. The K* mass turns 
out to be slightly heavier than experiment, though the 
supplemental [ir, K, <^]-input gives consistent results with 
experiment with large statistical error. Possible origin 
of the deviation is due to uncertainty of chiral fits. In 
fact, an alternative fit based on chiral perturbation the- 
ory (xPT) we discuss later yields wig, = 894(12) MeV 
(K- input). In Figj2] we overlay the previous results of 
meson masses Q in the Nf — 2 and quenched (Nf = 0) 
QCD with tadpole improved one- loop csw- The dynam- 
ical u and d quarks significantly reduce the O(10%) de- 
viation of the quenched spectrum from experiment. We 
find no further dynamical s quark effect beyond statisti- 
cal errors. 

The quark masses are evaluated for the MS scheme 
at the scale fi = 2GeV using the tadpole improved one- 
loop matching [l2T | at /i = a^ 1 with an improved cou- 
pling determined from plaquette and rectangular loop 
and four-loop renormalization group equation. In the 
continuum extrapolation of the quark masses, we assume 
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FIG. 3: Continuum extrapolations of the up, down and 
strange quark masses. For comparison, results for Nf = 
and N f =2 QCD @] are overlaid. 

the 0(g A ma) contributions are small and neglect it. As 
Fig. [3] shows, the quark masses are well described by 
a linear function in a 2 , and the values determined for 
either the K- or the ^-inputs, while different at finite 
lattice spacings, extrapolate to a common value in the 
continuum limit. Therefore the continuum limit is esti- 
mated from a combined linear fit with the K- and the 
0-inpute^ We obtain rn MS (fi = 2GeV) = 3.55(19) MeV 
and mf s {n = 2GcV) = 90.1(4.3) MeV with a suffi- 
ciently small x 2 /d.o.f. < 0.42. Note that the supplemen- 
tal [ir, K, </>]-input gives larger statistical error and hence 
is not used to estimate central values. 

We now turn to estimates of possible systematic errors. 
Finite size effect (FSE) — The meson masses at 
the infinite volume are estimated at /? = 1.90 using 
data on a V ~ (2.0fm) 3 lattice and those from our 
exploratory study on a V ~ (1.6fm) 3 lattice [l3j], and 
assuming a strong volume dependence of (mhad,v — 
)/mhad,v=oo oc I IV [14|. The chiral fits to the 
infinite volume values lead to less than a 4% change for 
the meson masses at the physical point. For the quark 
masses, however, we find a larger shift of 12.2% from a 
V ~ (2.0fm) 3 lattice to V = oo for m with <^-input and 
8.1% for m s with if -input (differences are smaller for the 
other cases). Assuming that FSE is independent of lat- 
tice spacing, we take the differences as estimates of FSE 
for the quark masses in the continuum limit. 
Chiral extrapolation — In addition to the polynomial 
chiral fits, we fit the meson masses using xPT formu- 
laes modified for the Wilson quark action (W%PT) [l5j ]. 
Namely, we fit m ff , m^., m 



p and m^, using the NLO 
Nf = 2 + 1 QCD WxPT formulae for the O(a) improved 
theory [lq |. Since the formula in Ref. [l6| is not appli- 
cable for the 4> meson, we estimate the effect only for 
if-input. In the fits we obtain m to be 3.1% smaller and 
m s to be 1.2% larger than those of the polynomial fit. We 
note that our W%PT fits to data do not exhibit a clear 
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chiral logarithm, probably because u and d quark masses 
in our simulation are not sufficiently small. Further pos- 
sible systematic error from a long chiral extrapolation for 
p , mentioned above, is estimated by the supplemental 
[tt, K, (/)]-input, which gives 3.0% larger for fn and 3.4% 
larger value for m s than the central one. For an estimate 
of systematic errors from chiral fits, differences of the two 
alternative fits from the central value are added linearly 
Renormalization factor — Uncertainty of the one- 
loop calculation of the renormalization factor is esti- 
mated by shifting the matching scale from p = 1/a to 
p = tt / a and also using an alternative tadpole improved 
coupling [2| . 

Continuum extrapolation — Possible 0(a 3 ) effects 
are investigated by performing the continuum extrapola- 
tion adding an 0(a 3 ) term to the fit function. 
Electromagnetic (EM) effects — Systematic error 
due to uncertainty of the EM effects is estimated fol- 
lowing extensive arguments 0, E3, EH to Dashen's 
theorem [Toj]. Namely, we estimate the effects by a 
further mass shift of our input using a relation 

(m 2 K± - m 2 K0 ) E M = (1 + A b )(to^ ± - to^ )exp assum- 
ing the EM effects for other mesons are negligeble. We 
vary the in range [—1,-1-1] as our estimate of the EM 
effects, and we find a quite small change in m s and no 
change in fn. 

Isospin breaking effects — Isospin breaking effects 
are estimated by chiral fits with Eq. ([I]) for m u ^ 
and taking to^o, m p o, m^± and m^-o as inputs. We 
find that m u /m^ = 0.577(25), and that m and m s 
have no change from the K input result. We note that 
fnu/fna strongly depends on an estimate of the EM ef- 
fects; m u /m d =0.663-0.498 for A E = [-1,+1], though 
fn and m s almost do not. 
Finally we obtain 

m m (n = 2GeV) 

= 3 1 55(19)(^ 3 )(±l 1 1 )(lf 7 )(^)(i°)(l°), (2) 
mf s (p = 2GeV) 

= 90.1(4.3)(^ 3 )(^ 2 )(^)(^ 2 - 8 )(^i)(^), (3) 

in MeV units, where the errors are statistical, systematic 
due to FSE, chiral extrapolation, renormalization factor, 
continuum extrapolation, EM effect and isospin break- 
ing effect, respectively. Adding the errors in quadrature 
yields the values quoted in the abstract. These values 
agree well with the latest report from the MILC Collab- 
oration |ra = 3.3± 0.3 MeV and m s = 90 ± 6 MeV 
where we added the quoted errors in quadrature. They 
also include the Nf = 2 values Q within the error. 

Scaling violation in the quark masses is unexpectedly 
large, while that for the meson masses are reasonably 
bounded at a percent level at a w 0.1 fm. To gain a 
better control over systematic uncertainties, a significant 
reduction in the simulated light quark masses on a 



correspondingly larger lattice is needed. An attempt is 
underway to meet these challenges (l9j . 
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